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ABSTRACT 



Context. The anelastic approximation is often adopted in numerical calculation with 



P4, 

, low Mach number, such as stellar internal convection. This approximation requires 

frequent global communication, because of an elliptic partial differential equation. 

Oh' Frequent global communication is negative factor for the parallel computing with a 

^ ■ large number of CPUs. 

5— i ' 

Aims. The main purpose of this paper is to test the validity of a method that artifi- 
cially reduces the speed of sound for the compressible fluid equations in the context 
of stellar internal convection. The reduction of speed of sound allows for larger time 
steps in spite of low Mach number, while the numerical scheme remains fully explicit 
and the mathematical system is hyperbolic and thus does not require frequent global 
communication. 

Methods. Two and three dimensional compressible hydrodynamic equations are solved 
numerically. Some statistical quantities of solutions computed with different effective 
Mach numbers (due to reduction of speed of sound) are compared to test the validity 
of our approach. 

Results. Numerical simulations with artificially reduced speed of sound are a valid ap- 
proach as long as the effective Mach number (based on the reduced speed of sound) 



$-H ' remains less than 0.7. 

Key words. SunJnterior - Sumdynamo - Method:numerical 



1. Introduction 

Turbulent thermal convection in the solar convection zone plays a key role for the mainte- 
nance of large scale flows (differential rotation, meridional flow) and solar magnetic activity. 
The angular momentum transport of convection maintains the global mean flows. Global 
flows relate to the generation of global magnetic field, i.e. the solar dynamo. Differential 
rotation bends the pre-existing poloidal field and generates the strong toroidal field (ft 
effect) and mean meridional flow transports the magnetic flux equatorward at the base 
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of the convection zone ( Choudhuri et al 



1995 



Dikpati fc Charbonneaul . [1999). The in- 



ternal structures of the solar differ ential rotation and the meridional flow are revealed by 



the helioseismology (see revie w by Thompso n et al 



2003). Some mean field studies have 



reproduced these global f l ow (IKichatinov fc Riidiger , 



1993: 



Ktiker & Stix 



2001 



2005 



Rempel . 



Hotta &: YokovamaL 1201111 . These studies, however, used some kinds of parameteri- 



zation of the turbulent convection, i.e. turbulent viscosity and turbulent angular momen- 
tum transport. Thus, a self-consistent thorough understanding of global structure requires 
the detailed investigation of the turbulent thermal convection. Turbulent convection is 
important also in the magnetic field itself. The strength of the next solar maximum in 
the prediction by the mean field model significantly depends on the turbulent diffusivity 
( Dikpati &: Gilmanl . 



2006 



Choudhuri et al 



2007 



Yeates et al 



20081 ). In addition, the tur- 



bulent diffusion has an important role in the pa rity of solar global field, the strength of 
polar field and so on ([Hotta fc Yokovama . 



2010 



ffl). 



T here are already numerous LES numerical simulations on the solar and stellar convec 



tion (Gilman 



Brown et al 



Brown et al 



1977 



2010 



Gilman fc Miller 



1981 



Glatzmaicr 




1984 






iesch et al. 




2000 


s (IGilman & Miller 




1981 




Brun et al. 



2006 



2004; 



20111 ). In these studies the anelastic approximation is adopted to avoid 



the difficulty which is caused by the high speed of sound. At the base of convection zone 
the speed of sound is ab out 200 km s _1 . In contrast the speed of convection is thought to 
be 50 m s _1 



(Stix 



2004), so the time step must be shortened due to the CFL condition in 



explicit fully compressible method, even when we are interested in phenomena related to 
convection. In the anelastic approximation, equation of continuity is treated as 

V • (pov) = 0, (1) 

where po is stratified background density and v denotes the velocity. The anelastic approxi- 
mation assumes that the speed of sound is essentially infinite, resulting in an instantaneous 
adjustment of pressure to flow changes. This is achieved by solving an elliptic equation for 
the pressure, which filters out the propagation of sound waves. As a result the time step 
is only limited by the much lower flow velocity. However, due to the existence of an el- 
liptic the anelastic approximation has a weak point. The numerical calculation in parallel 
computing requires the frequent global communication. At the present time the efficiency 
of scaling in parallel computing is saturated with about 2000-3000 CPUs in solar global 
simulation with pseudo-spectral method (M. Miesch private communication). More and 
more resolution, however, is thought to be needed to understand the precise mechanism of 
the angular momentum and energy transport by the turbulent convection and the behavior 
of magnetic field especially in thin magnetic flux tube. 

In this paper, we test the validity of a different approach to circumvent the severe 
numerical time step constraints in low Mach number flows. We use a method in which the 
spee d of sound is r educed artificially by transforming the equation of continuity to (see, 



Rempel 



2005) 



e.g. 

_ = --V>v), 



(2) 



where t denotes the time. Using this equation, the effective speed of sound becomes £ 
times smaller, but otherwise the dispersion relationship for sound wave remains unchanged 
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(wave speed is dropped for all wavelength equally). Since this technique does not change 
the hyperbolic character of the underlying equations, the numerical treatment can remain 
fully explicit and thus does not require global communication in parallel computing. We 
will call the technique in the following the Reduced Speed of Sound Technique (RSST). 
This technique has been used previously by Rempel (2005, 2006) in mean field models for 
solar differential rotation and non-kinematic dynamos, which essentially solve the full set 
of time dependent axisymmetric MHD equations. Those solutions were however restricted 
to the relaxation toward a stationary state or very slowly varying problems on the time 
scale of the solar cycle. Here we will apply this approach to thermal convection, where 
the intrinsic time scales are substantially shorter. To this end we study two and three 
dimensional convection, in particular the latter will be non-stationary and turbulent. 

The detailed setting of test calculation is given in Section The results of our calcula- 
tions are given in Section [3] We summarize our paper and give discussion of the RSST in 
Section g] 



2. Model 

2.1. Equations 

The two or three dimensional equation of continuity, equation of motion, equation of energy, 
equation of state, are solved in Cartesian coordinate (x, z) or (x,y,z), where x and y 
denote the horizontal direction and z denotes the vertical direction. The basic assumptions 
underlying this study are as follows. 

1. Time independent hydrostatic reference state. 

2. The perturbations caused by thermal convection are small, i.e. p\ <C po and p\ <C po, 
Here po and Po denote the reference state values, whereas p\ and p\ are the fluctuations 
of density and pressure, respectively. Thus a linearized equation of state is used as eq. 
@ 

3. The profile of the reference entropy sq(z) is a steady state solution of the thermal 
diffusion equation V • {K pqTqV sq) — with constant K . 



The formulations are almost same as 



Fan et al 



(|2003l ). Equations are expressed as, 



^ = -^V.( Po v), (3) 

| = -(v.v)v-^-^ + iv.n, (4) 

ot po po po 

d 1 

= -(v ■ V)(« + «i) + -^V • (Kp ToV Sl ) 
ot poTq 

' ' II V: v. (5) 



Po 



Pi =Po [7— + fi i J i ( 6 ) 
V Po / 

where Tq(z), and sq(z) denote reference temperature, and entropy, respectively and e z 
denotes the unit vector along the z-direction. 7 is the ratio of specific heats, with the 
value for an ideal gas being 7 = 5/3. s± denotes the fluctuation of entropy from reference 
atmosphere. Note that the entropy is normalized by specific heat capacity at constant 
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volume c v . The quantity g is the gravitational acceleration, which is assumed to be constant. 
The quantity II denotes the viscous stress tensor, 



Iljj = p v 



dx-j 



dvj_ 

dxi 



(?) 



and v and K denote the viscosity and thermal diffusivity, respectively, v and K are assumed 
to be constant throughout the simulation domain. We assume for the reference atmosphere 
a weakly superadiabatically stratified polytrope: 

z 



Po(z) = p. 
Po{z 



1 



Pr 



(to + l)H r 

z 



-i m+l 



(TO + f )H r 



T Q (z) = T r 



1 - 



(m + f )H r 



H p (z) = 

ds 
dz 

5{z) 



Pa 

Pr 

r po(z)' 



(8) 

(9) 
(10) 
(11) 
(12) 
(13) 



where p t ,Pr, T T , H T , and 5 T denote the values of po, po, Tq, Hq (the pressure scale height) and 
5 (the non-dimensional superadiabaticity) at the bottom boundary z — 0. Since |<5| <C 1, 
the value m is nearly equal to the adiabatic value, meaning to = 1/(7 — 1). The strength of 
the diffusive parameter v and K is expressed with following non-dimensional parameters: 
the Reynolds number Re = v c H T /i>, and the Prandtl number Pr = v / K , where the velocity 
unit i> c = (^S^gHf) 1 / 2 . Note that in this paper the unit of time is H r /v c . In all calculations, 
we set Pr = 1. 



2.2. Boundary Conditions and Numerical Method 

We solve equations ©-© numerically. At the horizontal boundaries (x — 0, L x and 
y = 0, Ly), periodic boundary conditions are adopted for all variables. At the top and the 
bottom boundaries impenetrative and stress free boundary conditions are adopted for the 
velocities and the entropy is fixed: 

v z = 0, (14) 

ar = ' (15) 

7T*=0, (16) 
oz 

Sl = 0. (17) 

At both top and bottom boundaries (z — and z — L z ), we set pi in the ghost cells 
such that the right hand side of the z-component of eq. (j4|) is zero at the boundary (which 
is between ghost cells and domain cells), where the ghost cells are the cells beyond the 
physical boundary. 

We adopt the fourth-order space-centered difference for each derivative. The first spatial 
derivatives of quantity q is given by 
fdq\ 1 



dx / , 12Ax 



(-<?i+2 + 8q i+ i - 8^-1 + qi-2), (18) 
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where i denotes the index of the grid position along a particular spatial direction. The 
numerical solution of the system is advanced in time with an explicit fourth-order Runge- 
Kutta scheme. The system of partial equations can be written as 
<9U 

^=R(U) (19) 



(20) 

(21) 

(22) 
(23) 

The maximum allowed time step, Ai max , is determined by the CFL criterion. When 
both advection and diffusion terms are included in calculation, the time step reads, 

At max = mm(At ad ,A* v ). (24) 

Here 

A . min(Ax, Ay, Az) 

At ad = C ad • (25) 

Ctot 

ctot is the total wave speed: 

ctot = \v\ + c' s , (26) 
where the effective speed of sound is expressed as: 



n+i, which is 


the value at t n +i 


U n+ i = u„ 


+ ^R(U n ), 


U„+i = u„ 


+ ^R(u„+i), 




+ ^R(U n+ i), 


U„+i = U„ - 


|-A*R(U n+ i). 



< = lhf- (27) 

The time step determined by diffusion term is 

_ min(Ax 2 ,Ay 2 ,Az 2 ) 
i±i v — c v . . yzaj 

max(iv, v) 

c a d and c v are safety factors of order unity. 

Using S T = 1 X 10" 4 , the original speed of sound is about 35w c at the bottom and 12v c 
at the top boundary, respectively. In all calculation Ax ~ 2.3 x 10~ 2 iJ r . Thus if we use 
c ad = 1 and c v = 1, At ad = 6.4 x 1CT 4 H T /v c and At v = 1.5 x 10 _1 H T /v c . For all the ^ 
values considered in this paper the time step remains restricted by the (reduced) speed of 
sound, thus the calculation has about £ times better efficiency with the RSST. 



3. Result 

3.1. Two Dimensional Study 

Four two-dimensional calculations with the RSST and one calculation without approxima- 
tion are carried out (case 1-5). The values of free-parameters are given in Table [TJ The 
superadiabaticity (5) is 1 x 10~ 6 at the bottom and about 2 x 10~ 5 at the top boundary. 
It is almost the same superadiabaticity value as the base of solar convection zone. Fig. Q] 
shows the time-development of entropy. In the beginning the non-linear time dependent 
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convection developed (top panel), which transitions to a steady state at later times (bot- 
tom panel). In the steady state, we compare the RMS velocity with different £. The RMS 
velocity is defined as: 



Fig- HI shows our results with the value of £ being from 1 to 80. The effective Mach number 
is defined as Ma = fRMs/Cg. If we have larger £ value than 80, we cannot obtain stationary 
state, since there are some shock generated by supersonic convection. The discussion of 
unsteady convection is given in the next session of three-dimensional calculations. Even 
though the Mach number reaches 0.6 using £ = 80 (panel a), the horizontal and vertical 
RMS velocity is almost same as those calculated with £ — 1 (without RSST). The ratio 
between the RMS velocities with each £ and £ = 1 are shown in Fig. [2^. The deviation is 
always a few percent. This result is not surprising since in the stationary state the equation 
of continuity becomes 



whose solution must not depend on the value of £. We note that the cell size is to some 
degree affected by the aspect ratio of the domain. This does not affect our conclusions since 
this influence is the same for all values of £ considered. To confirm the robustness of our 
conclusions we repeated this experiment with a wider domain (26.16iJ r x 2.1SH r instead 
of 8.72H r x 2. 18-ff,-) in which we find 4 steady convection cells with about 10% different 
RMS velocities. Also here we find a dependence on £ similar to that shown in Fig. [5J i.e. 
the solutions show differences of only a few percent as long as £ < 80. In this section we 
confirm that the RSST is valid for the two-dimensional stationary convection when the 
effective Mach number is less than 1, i.e. £ < 80 with S T = 1 X 10~ 6 . If we use 6 T = 1 x 10~ 4 
(result is not shown), the criterion becomes £ < 8 in two dimensional calculation. 

As a next step we investigate the dependence of the linear growth rate on £ during the 
initial (time dependent) relaxation phase toward the final stationary state. Fig. [3] shows 
the linear growth of maximum perturbation density p\ with different £. Black and red lines 
show the results with 5 t = 1 x 10~ 6 and 6 T = 1 x 10~ 4 , respectively. These calculation 
parameters are not given in Table Q] In the calculation with 8 = 1 X 10~ 6 the growth 
rate decreases for values of £ > 300, whereas it occurs £ > 30 in the calculation with 
6 T = 1 X 10~ 4 . The reason can be explained as follows: In the convective instability, upflow 
(downflow) generates positive (negative) entropy perturbation and then negative (positive) 
density perturbation is generated by the sound wave. If the speed of sound is fairly slow, 
the generation mechanism of density perturbation is ineffective. In the calculation with 
5 T = 1 x 10~ 6 (<5 r = 1 x 10~ 4 ), the growth rate with £ = 200 (£ = 20), however, is almost 
same as that with £ = 1, even though flow with £ = 200 (£ = 20) is expected to be the 
Mach number Ma = 1.3, i.e. supersonic convection flow in the saturated state. 

3.2. Three Dimensional Study 

In this section, we investigate the validity of the RSST with three-dimensional unsteady 
thermal convections (case 6-13). The value of superadiabaticity at the bottom boundary is 




(29) 



= 72 V • fov), 



(30) 
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1 x 1CP 4 and at the top 2 x 10~ 3 . Although this value is relatively large compared with solar 
value, the expected speed of convection is much smaller than speed of sound, so it is small 
enough to investigate the validity of the RSST. Entropy of three-dimensional convections 
with £ = 1, 20, and 80 at t = \QQH r /v c are shown in Fig. SJ The convection is completely 
unsteady and turbulent (animation is provided). The appearance of convection with £ = 20 
is almost same as that with £ = 1. This will be verified below by the Fourier transformation 
and auto-detection technique. On the other hand, the appearance of convection with £ = 80 
(bottom panel) is completely different from the others. This difference is best visible in the 
animation of Fig. 0] that is provided with the online version. 

RMS velocities with different £ are estimated as average of values between t = 100 to 
200ff r /?; c (Fig. [5J panel b and c). Without the RSST, i.e. £ = 1 and using 5 r = 1 x 10~ 4 , 
the Mach number is 1 x 10~ 2 at the bottom and 4 x 10~ 2 at the top boundary. The 
RMS velocities at £ = 40 and 80 differ from those at the £ = 1 by less than 15 % and 
30 %, respectively. When we adopt £ = 40 and 80, the Mach number estimated by RMS 
velocity exceeds unity (Fig.[SJ panel a), i.e. supersonic convection. This supersonic downflow 
frequently generates shocks and positive entropy perturbation, thus downflow is slowed. 
This is the reason why the RMS velocities with large £(= 40, 80) are small. When £ = 5, 
10, 15, and 20, however, the RMS velocities show good agreement with that with £ = 1. 
RMS power density of pressure, buoyancy and inertia are estimated (Fig. [5]). This results 
shows almost the same tendency as the result of RMS velocities. Power profiles £ = 5, 
10, 15, 20 show agreement and those with £ = 40, 80 shows discrepancy with those with 
£ = 1 . These results show that the RSST is valid technique with at least £ = 20 with which 
the Mach number is around 0.7. Note that the convection pattern in our 2D and 3D cases 
differs substantially. As a consequence, also the £ values for which the validity of RSST 
breaks down are also different in 2D and 3D setups. 

In Fig. [6] we compare averaged spectral amplitudes for different values of £. There is no 
significant difference between the spectral amplitudes for different values of £ in the range 
from 1 to 20. 

We investigate the distribution of cell size at the top boundary. This value is significantly 
related to the turbulent diffusivity and transport of angular momentum or energy. The 
method to detect the cell is explained as follows. At the boundary of the cell i.e., the 
region of downflow, the perturbation of density is positive and has large value. When the 
density in a region exceeds a threshold, the region is regarded as boundary of convective 
cell. When a region is surrounded by one continuous boundary region, the region is defined 
as one convective cell. Fig. [7] shows the detected cells. Each color and each label (#n) 
correspond to each detected cell. We estimate size of all cells and compare the distribution 
of cell size with different £. The results are shown in Fig.[5J The cell size distribution follows 
a power law from about 0.01 to 10 and there is no dependence on £ in the range from 1 
to 20. Although using this type of technique size of cells is tend to be large with neglecting 
smaller cells, our conclusion is not wrong, since all the auto detections are affected equally. 

With above two investigations, i.e. the Fourier analysis and study of cell detection, we 
can conclude that the statistical features are not influenced by the RSST as long as the 
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effective Mach number (computed with the reduced speed of sound) does not exceed values 
of about 0.7, which corresponds to £ = 20 in our setup. 

In order to confirm our criterion that the RSST is valid if the Mach number is smaller 
than 0.7, we conduct calculations with larger superadiabaticity, i.e. <5 r = 1 x 10~ 3 (Case 
14-18). If our criterion is valid, required £ with larger superadiabaticity must decrease. The 
results are shown in Fig. O Using 6 t = 1 X 10 -3 , the calculations with £ = 10, 15, and 20 
generate supersonic convection flow near the surface (Fig. |Hk). It is clear that the results 
with £ — 10, 15 and 20 differ from that with £ = 1 and 5. The calculation with £ = 5 
shows flow whose Mach number is around 0.6. Thus our criterion is not violated even with 
different value of the superadiabaticity. 

Due to our changed equation of continuity the primitive and conservative formulation 
of Eq. ([3]) to (J5J) are not equivalent anymore. For example, the equation of motion in 
conservative form is expressed as: 

^(pv) = -V-(pvv)+F, (31) 

where F denotes pressure gradient, gravity and Lorentz force. With some transformations 
we can obtain, 

v^+^ = -vV-(pv)-p(vV)v + F (32) 

If equation of continuity is satisfied, i.e. {dp/dt = — V • (pv)), the primitive form is obtained 
as: 

p— = p(v-V)v + F. (33) 

However, using the RSST these two form are no longer equivalent. We used here the 
primitive formulation at the expense that energy and momentum are not strictly conserved; 
however, the consistency of our results for different values £ strongly indicates that this is 
not a serious problem for the setup we considered. Alternatively we could also implement 
our modified equation of continuity into a conservative formulation. This would ensure that 
density, momentum and energy are strictly conserved at the expense of a modified set of 
primitive equations. Fig. [TU] shows the dependence of [V • (pov)]rms (= [£ 2 <9pi/<9i]RMs) 
on £. Using £ = 5, this term is almost same as that with £ = 1. Although the deviation 
becomes large as £ increases, it is not proportional to £ 2 . 

In the previous discussion we kept £ constant in the entire computational domain. In 
the solar convection zone the Mach number varies however substantially with depth, from 
~ 1 in the photosphere to < 10~ 7 at the base of the convection zone, the speed of sound 
itself varies from about 7 km s _1 in the photosphere to 200 km s _1 at the base of the 
convection zone. A reduction of the speed of sound is therefore most important in the deep 
convection zone, but not in the near surface layers. This could be achieved with a depth 
dependent £. Even if we use conservative form as, 

•(?")• <*> 

the result must not be same as the result without approximation in statistical steady state. 
When the value averaged in statistical steady state is expressed as (a) , where a is a physical 
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= V-(-p (v>), (35) 



value, the equation of continuity becomes 
1 

I 2 

with inhomogeneous £. The solution of Eq. (|35[) is different from the solution of original 
equation of continuity in statistical steady state, i.e. = V • (po(v)). Thus, the statistical 
features such as RMS velocity are not reproduced with inhomogeneous £. If we use the 
non-conservative form of equation of motion as Eq. @, this problem does not occur. Thus 
we investigate the nonuniformity of £ in case 13 using non-conservative form, i.e., eq. ©. 
In order to keep the Mach number uniform in all the height, we use £ = 20/(<5(z)/(5 r ) 1 / 2 , 
i.e. £ = 20 at the bottom and £ = 4.5 at the top boundary. Note that the ratio of the speed 
of convection to the speed of sound is roughly estimated as y/S. In this setting the value 
of y/S is 1 x 10~ 2 at the bottom and 4 x 10 -2 at the top boundary. The same analysis 
as that for uniform £ is done (see Fig. [5J [5] and ITU)) . There is no significant difference 
between £ = 20/ (5(z)/5 r ) 1 ' 2 and £ = 1. Although mass is not conserved locally with 
non-homogeneous £ using non-conservative form, horizontally averaged vertical mass flux 
is approximately zero in the statistically steady convection and the conservation total 
mass is not significantly broken. We conclude that an inhomogeneous £ is valid under the 
previously obtained condition, i.e. the Mach number is less than 0.7. 



4. Summary and Discussion 

In this paper we applied RSST (see Section [TJ to two and three dimensional simulations of 
low Mach number thermal convection and confirmed the validity of this approach as long 
as the effective Mach number (computed with the reduced speed of sound) stays below 0.7 
everywhere in the domain. The overall gain in computing efficiency that can be achieved 
depends therefore on the maximum Mach number that was present in the setup of the 
problem. 



Since the Mach number is estimated to be 10 4 in the base of solar convection zone ( Stix 



2004), several thousand times longer time step can be taken using the RSST. Therefore the 
RSST and parallel computing with large number of CPUs will make it possible to calculate 
large scale solar convection with high resolution in the near future. 

Compared to the anclastic approximation there are three major advantages in RSST: 

1. It can be easily implemented into any fully compressible code (regardless of numerical 
scheme or grid structure) since it only requires a minor change of the equation of 
continuity and leaves the hyperbolic structure of the equations unchanged. 

2. Due to the explicit treatment it does not require any additional communication over- 
head, which makes it suitable for massive parallel computations. 

3. The base of the convection zone and the surface of the sun where the anelastic approx- 
imation is broken can be connected, using space dependent £. 

Overall we find that RSST is a very useful technique for studying low Mach number flows in 
stellar convection zones as it substantially alleviates stringent time step constraints without 
adding computational overhead. 
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Table 1. Parameters for numerical simulations. In case 13 the value of £ is 20 at the bottom 
and 4.5 at the top boudary. 



Case 


dimension 


L x x L 


y X Z/ 2 (_ff r ) 


N x x N y : 


x N z 




S r 




Re 


£ 


1 


2 


8 


.72 x 2.18 


384 x 96 


1 


x 10" 


-6 


260 


l 


2 


2 


8 


.72 x 2.18 


384 x 96 


1 


x 10" 


-6 


260 


10 


3 


2 


8 


.72 x 2.18 


384 x 96 


1 


x 10" 




260 


30 


4 


2 


8 


.72 x 2.18 


384 x 96 


1 


x 10" 




260 


50 


5 


2 


8.72 x 2.18 


384 x 96 


1 


x 10" 


-6 


260 


80 


6 


3 


8.72 


X 


8.72 


X 


2.18 


384 


X 


384 


x 96 


1 


x 10" 




300 


1 


7 


3 


8.72 


X 


8.72 


X 


2.18 


384 


X 


384 


x 96 


1 


x 10" 




300 


5 


8 


3 


8.72 


x 


8.72 


X 


2.18 


384 


X 


384 


x 96 


1 


x 10" 




300 


10 


9 


3 


8.72 


x 


8.72 


X 


2.18 


384 


X 


384 


x 96 


1 


x 10" 




300 


15 


10 


3 


8.72 


x 


8.72 


X 


2.18 


384 


X 


384 


x 96 


1 


x 10" 




300 


20 


11 


3 


8.72 


x 


8.72 


X 


2.18 


384 


X 


384 


x 96 


1 


x 10" 




300 


40 


12 


3 


8.72 


x 


8.72 


X 


2.18 


384 


X 


384 


x 96 


1 


x 10" 


-4 


300 


80 


13 


3 


8.72 


X 


8.72 


X 


2.18 


384 


X 


384 


x 96 


1 


x 10" 


-4 


300 


20/{5/8 r ) 1/2 


14 


3 


8.72 


X 


8.72 


X 


2.18 


384 


X 


384 


x 96 


1 


x 10" 


-3 


300 


1 


15 


3 


8.72 


X 


8.72 


X 


2.18 


384 


X 


384 


x 96 


1 


x 10" 


-3 


300 


5 


16 


3 


8.72 


X 


8.72 


X 


2.18 


384 


X 


384 


x 96 


1 


x 10" 


-3 


300 


10 


17 


3 


8.72 


X 


8.72 


X 


2.18 


384 


X 


384 


x 96 


1 


x 10" 
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Entropy (8 r c v ) -^ = 0.50 




Fig. 1. Time-development of entropy in a two-dimensional calculation with parameters of 
case 1. (Top panel) t = 20. (Bottom panel) t = 400. 



(a) Maximum Mach number (b) Mach number 




z (H r ) z (H r ) z (H r ) 

Fig. 2. Some quantities in two dimensional calculation. (a) Maximum Mach number of 
each time step, (b) Distribution of Mach number estimated with the RMS velocity, (c) 
Distribution of horizontal RMS velocity «h- (d) Distribution of vertical RMS velocity, (e) 
Ratio between the RMS velocities with each £ and £ = 1. 
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Fig. 4. Snapshots of entropy of the three-dimensional convection. Top, middle, and bottom 
panels correspond to £ = 1, 10 and 80 respectively. Left (right) panels show entropy at top 
(bottom) boundary. (Animation is provided, and the difference between £ — 1 and 80 is 
best visible in the animation of Fig. 0] that is provided with the online version.) 
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(a) Mach number (b) Vertical RMS velocity 
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(c) Horizontal RMS velocity 
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Fig. 5. Some quantities averaged in time between t = 100 and 200H r /v c in three dimen- 
sional calculation, (a) Distribution of Mach number estimated with the RMS velocity, (b) 
Distribution of vertical RMS velocity Vh- (c) Distribution of horizontal RMS velocity, (d) 
Distribution of RMS power density of pressure and buoyancy, (e) Distribution of RMS 
power density of inertia. 
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Fig. 6. Comparison of horizontal velocity spectra with different £. (a) z = 2.18H r (b) 
z = l.09H r (c) z = 0. Velocities are averaged in time between t — 100 and 200H r /v c 
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(a) at z=2.1 8H r (b) Detection of convective cell 
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Fig. 7. Detection of convective cell, (a) Original contour of density perturbation at 
z = 2.18H r (b) Distribution of detected convection cell. Color and label (#n) show each 
convective cell. 



Cell Size Distribution 
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Fig. 8. Distribution of convective cell size with different £ is shown. The cell size distribu- 
tion is averaged in time between t = 100 and 200H r /v c . 
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Fig. 9. The results with S T = 1 x 10 3 . The format is the same as Fig. [5] 
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Fig. 10. Dependence of [V • (po v )]rms on £. In case 13, the value of £ is 20 at the bottom 
and 4.5 at the top boundary. 
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